***********************************************************************************************************
* This replication file maps the notary locations and municipal population density to replicate Figure A1 *
***********************************************************************************************************
* 1. Save the locations on the notaries in a format which is appropriate for mapping the information.
* 2. Save the administrative data on census unit boundaries in a format which is appropriate for mapping the information.
* 3. Merge the GIS data on census units with population statistics to calculate population density.
* 4. Plot the notary locations and municipal population density to replicate Figure A1.
*************************
* 0 Set path to files *
*************************
clear all
capture log close
set more off
set trace off
cd $main_directory


*******************************************************************************************************
* 1. Save the locations on the notaries in a format which is appropriate for mapping the information. *
*******************************************************************************************************

use "Data_Belgium\Belfirst_data.dta", clear
keep off_id Lat_ext Lon_ext
rename Lon_ext _X
rename Lat_ext _Y
local lon0 "4.476940500000069"
local f "298.257223563"
local a "6378137"
geo2xy _Y _X, replace proj(mercator, `a' `f' `lon0')
save "Data_Belgium\GIS_data_census_units\geo2xy_notary_coordinates.dta", replace

****************************************************************************************************************************
* 2. Save the administrative data on  census unit boundaries in a format which is appropriate for mapping the information. *
****************************************************************************************************************************

clear
use "Data_Belgium\GIS_data_census_units\units_coord.dta"
rename _ID id
merge m:1 id using "Data_Belgium\GIS_data_census_units\units_db.dta"
gen country=substr(CENSUS_ID,1,2)
//Keep only the observations in Belgium
keep if country=="BE"
rename id _ID
keep _ID _Y _X
sort _ID
geo2xy _Y _X, replace proj(mercator)
save "Data_Belgium\GIS_data_census_units\geo2xy_census_coor_proj.dta", replace

*****************************************************************************************************
* 3. Merge the GIS data on census units with population statistics to calculate population density. *
*****************************************************************************************************

use "Data_Belgium\GIS_data_census_units\units_db.dta", clear
//Generate a country identifier
gen country=substr(CENSUS_ID,1,2)
//Keep only the observations in Belgium
keep if country=="BE"
sort id
//Generate the NIS variable
gen nis=substr(CENSUS_ID,7,5)
destring nis, replace
merge 1:1 nis using "Data_Belgium\Generated_data\BE_municipal_data.dta", force
drop _merge
gen population_density=population/Area_sqkm

*****************************************************************************************
* 4. Plot the notary locations and municipal population density to replicate Figure A1. *
*****************************************************************************************

colorpalette Greys, n(7) nograph
local colors `r(p)'
//Figure A1
spmap population_density using Data_Belgium\GIS_data_census_units\geo2xy_census_coor_proj, id(id) fcolor("`colors'") clmethod(custom) clbreaks(0 250 500 1000 2000 30000) legend(size(medium)) polygon(data("Data_Belgium\GIS_data_census_units\geo2xy_census_coor_proj") ocolor(black) osize(0.08) legenda(off)) point(data("Data_Belgium\GIS_data_census_units\geo2xy_notary_coordinates") x(_X) y(_Y) shape(D ..) fcolor(black) ocolor(black) osize(thin) size(0.5 ..) legenda(on) leglabel("notary offices")) legend(label(6 "2000+"))
graph save Estimates\Graphs\figure_A1, replace
graph export "Estimates\Graphs\figure_A1.eps", replace
